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We study a quasi-2D classical Landau-Ginzburg-Wilson effective field theory in the presence of 
quenched disorder in which incommensurate charge-density wave and superconducting orders are 
intertwined. The disorder precludes long-range charge-density wave order, but not superconducting 
or nematic order. We select three representative sets of input parameters and compute the corre¬ 
sponding charge-density wave structure factors using both large-A techniques and classical Monte 
Carlo simulations. Where nematicity and superconductivity coexist at low temperature, the peak 
height of the charge-density wave structure factor decreases monotonically as a function of increas¬ 
ing temperature, unlike what is seen in X-ray experiments on YBa2Cu306-i-a;. Conversely, where the 
thermal evolution of the charge-density wave structure factor qualitatively agrees with experiments, 
the nematic correlation length, computed to one-loop order, is shorter than the charge-density wave 
correlation length. 


I. INTRODUCTION 

The cuprate superconductors manifest remarkably 
rich phase diagrams, but with many features common 
to different materials within the familyIn ad¬ 
dition to the well known antiferromagnetic insulating 
and superconducting (SC) phases, recent experiments 
have revealed that short-range-correlated incommensu¬ 
rate charge-density wave (CDW) order, long known to 
play a prominent role in the physics of a limited sub¬ 
set of cuprates^^^, occurs in one way or another, in all 
or most cuprates®^^®. The thermal evolution of the X- 
ray structure factor^®“^®d 7 , is, 20-22 these newly stud¬ 

ied cases typically exhibits a gentle “concave-upward” 
onset rather than the sharp onset that is expected at 
the point of a thermodynamic phase transition. Mean¬ 
while, transport^®^^®, neutron scattering^®, 
and NMR®^ measurements have revealed anisotropies 
that are suggestive of the existence of long-range ne¬ 
maticity (broken C 4 rotation symmetry) in an overlap¬ 
ping regime of the phase diagram. 

The present theoretical study is carried out with ob¬ 
servations in the model cuprate YBa 2 Cu 306 +a; (YBCO) 
in mind, specifically for a range of temperatures, T, low 
compared to the “pseudo-gap” crossover®’®, T*, and for a 
range of doping concentrations 5 in the neighborhood of 
i5 = 1/8 where CDW fluctuations are experimentally de¬ 
tectable. In YBa 2 Cu 3 06.67 (corresponding to i5 « 0.12), 
the observations of Ref. 14 then restrict us to T < 150 K. 
We do not comment here on the higher T regime where 
there are thermodynamic or spectroscopic indications of 
a pseudogap. We will use a fluctuating order model to 
study the T dependence of the intertwined CDW and SC 
orders, and their relation to nematicity. 

Previous theoretical works®®^"'® introduced classical 
Landau-Ginzburg models to study the competition 
among different order parameters in cuprates. In par¬ 
ticular, Ref. 35 focused on the angular fluctuations of a 


multi-component order parameter, consisting of SC and 
CDW correlations in two spatial dimensions without dis¬ 
order. Ref. 36 investigated the effects of quenched dis¬ 
order and dimensionality on CDW and nematic orders 
in the cuprates. Here, we consider a generic Landau- 
Ginzburg theory with a multi-component order parame¬ 
ter (consisting of one SC complex field 4/ and two CDW 
complex fields ^x,y) and quenched disorder in a quasi-2D 
system. We show that the T-dependence of the CDW 
structure factor depends strongly on the strength of the 
disorder, the dimensionality of the system, and also on 
other input parameters of the model. We also calculate 
other quantities such as the nematic correlation length 
and the integrated intensities of SC and CDW orders. 

While our work was being completed, we learned of 
the similar analysis by Caplan et They tune the 

competition between CDW and SC by an applied mag¬ 
netic held, and obtain trends consistent with our results 
in Sec. IV A below. 

The format of this paper is as follows. In Sec. H, we 
introduce a classical Landau-Ginzburg model for a lay¬ 
ered system with SC and CDW orders and random-held 
type disorder. Sec. HI illustrates the methods we use to 
solve this model, including the replica trick and a large- 
N expansion, which are applied to obtain a saddle-point 
(mean-held) solution of the model, as well as classical 
Monte Carlo techniques. In Sec. IV, we report results for 
the T-dependence of the CDW structure factor in various 
regions of the phase diagram, and we discuss the effects of 
both disorder and dimensionality. We present in Sec. V 
detailed mean-held phase diagrams as functions of vari¬ 
ous input parameters and temperature, and in Sec. VI, 
we show calculations of the nematic correlation length 
to one-loop order in the non-nematic phase. Finally, in 
Sec. VH, we discuss the implications of our results as well 
as connections to cuprate experiments. 
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II. THE MODEL 


We consider a layered system with tetragonal symme¬ 
try. The charge density at position r in layer m can be 
expressed as 


Pm(i’) = P + -h c.c.] 


( 2 . 1 ) 


where p is the uniform charge density, Qx and Qy are 
incommensurate in-plane wave vectors along the x and 
y directions with equal magnitude, and $ 3 , and are 
classical CDW order parameters with slow spatial vari¬ 
ation. The model we study is an effective field theory 
of $ = and an SC order parameter in the 

presence of random-field disorder h. The corresponding 
classical Hamiltonian is 


H = - +C.C. 

(r,r') ,m 

- + ^44'L(r)^'m-n(r) -f c.c. 

r,m 

- x) - $LW'r^’m(r+ y) + C.C. 

r,?7i 

r,m 

r,?7i 

+ E 

r,m 

+ E [^m(i’)^m(r) -f C.C.] , (2.2) 


where (dagger here denotes conjugate trans¬ 

pose of a complex vector), the lattice constant is set equal 
to 1 , r and r' are nearest-neighbor xy-plane coordinates, 
X and y are ccy-plane unit vectors, m labels the layer 
along z direction, N is the number of real components of 
each order parameter, and 


relates the six real components of (d>, $ 2 ,, $j,). More gen¬ 
erally, in the absence of the random field, the model has 
a symmetry of SO{2) x S'0(2) x SO[2) x Z 2 , where 
comes from invariance under 7 r /2 spatial rotation. The 
random field breaks the symmetry to a single SO{2) (cor¬ 
responding to the phase of the SC order parameter), al¬ 
though the remaining symmetries are respected on aver¬ 
age in the disorder ensemble. 

Thermodynamic stability requires C/ > 0. To simplify 
the analysis, we consider the limit U — )■ -foo, which is 
equivalent to imposing the constraint 

|$x|" + |$yP + |^P = 3fV. (2.4) 

This constraint is a reflection of the experimental evi¬ 
dence that SC and CDW compete with each other at low 
temperature^®“^®d 7 , 20 , 22 _ Hamiltonian (2.2) then be¬ 
comes a non-linear sigma model. Note, however, that cal¬ 
culations can be carried out in the same manner for large 
but finite U, which give qualitatively similar results. 

The constraint in Eq. (2.4) leads to an equivalency be¬ 
tween Eq. (2.2) and the Hamiltonian studied by Monte 
Carlo methods in Refs. 35 and 37. Appendix A pro¬ 
vides information about how to relate parameter values 
used in the present and previously studied models. Note, 
however, that these previous works did not consider the 
effects of random-field disorder h and interlayer couplings 
Jz and 14 . Such effects will be studied in detail in this 
paper. (These previous studies also excluded from the 
Hamiltonian the term proportional to J', but the effects 
of this term were discussed in detail in Ref. 36.) 

g and g' determine the relative energy cost of order¬ 
ing between CDW and SC. The sign of A distinguishes 
between stripe (unidirectional CDW) and checkerboard 
phases. In our calculation we always take A > 0 (fa¬ 
voring stripes). Note that, although it does not break 
C 4 rotational symmetry, a positive A favors anisotropy 
between CDWs along x and y directions. 

The disorder potential h is taken to be a Gaussian 
random field with 


— Oj 


(2.5) 


IN/ 2 .N /2 \ (2.3) 

— ^N/2xN/2 J 

under the basis -|- ,..., -|- -f 

..., -I- . All of the following Monte 

Carlo results set N = 2 (as in Refs. 35 and 37). The terms 
proportional to J, AT, and Vz are the lattice versions 
of the familiar gradient terms in the continuum. Our 
large-fV and Monte Carlo calculations take J = K and 
Jz = Vz- (See Sec. VHA for discussion of Jz yf 14■) 
Because we are interested in quasi-2D systems, we always 
consider the case where 1 >> J^/J > 0. 

For the special case in which J — K and Jz = 14, 
with g = g' = J' = A = 0, and in the absence of a ran¬ 
dom field, the model has a large SO(6) symmetry which 


and 

I* \ (^■^) 

where denotes a disorder configuration average, 
a,/3 = x,y and i,j = 1,...,A. Notice that any lin¬ 
ear couplings between h and the SC order parameter dt 
are forbidden by gauge invariance. 

III. METHODS 

A. Saddle-point solution in the large-A limit 

We apply the replica trick^^ and integrate out h, 
then decouple the quartic terms using two Hubbard- 















consistently are 
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Stratonovich (HS) auxiliary fields. The resulting Hamil¬ 
tonian is 


.^replica 


- +C.C. 

(r.r'> 

m,a 

- E [j.‘J>E(r)$a,™+i(r) + y. m(^)^a,m+l (r) + C.C. 
r,m,a 

r,m,a 

+ E *Aia,m(r) |^'a.m(r)P + |^'a.m(r)P - 3A^ 

r,m,a 

- E (5 + 6 ff')l^a,m(r)P 

r,m,a 

+ E [»^a,m(r)|^a.m(r)P + 

r,m,a ^ 

T ^ ^ 771(1*) I ^a,a:,777(1*) I | d^a,y, 7 T 7 (l*) | 


r,r77,a 

20 ^ 

T 


■A'a,m(l*)Af 

4A 


E ^ll.m(r)$a 2 .m(r) 


(3.1) 


r,m 
ai ,02 


where 0 , 01,02 are replica indices; the HS fields M and 
rj in the second and third to last lines correspond to the 
quartic terms proportional to A and g' in Eq. (2.2) re¬ 
spectively; A/' 77 i(r) is the nematic order parameter. The 
Lagrange multiplier g enforces the constraint in Eq. (2.4). 
The expectation values (averaged both thermally and 
over disorder realizations) of CDW and SC order param¬ 
eters are obtained via diagonalization of Eq. (3.1) in both 
Eourier and replica spaces. Assuming there is no replica 
symmetry breaking^^, we obtain 

(riW.l.(r»^iv/4A(£_ + |!),,3.2, 
(4.l(r)4.,(r)) = iV• <«) 

_ r 7 / 3 u 77 

(^t(r)vk(r)) = nJ —(3.4) 

where (...) and denote thermal and disorder averages 
respectively, and 

= —2J{coskj; + cosky)^2J'{coskj: — cosky) 

— 2 cos fcz ± Af-I-(3.5) 

C = —2K{cos kx + cos ky) — 2Vz cos kz 

- g-6g' + g + g, (3.6) 

where we have taken the mean-field approximation that 
Aim(i*) = = Af, ??777 (r) = V for any (m,r). 

The saddle-point equations that are to be solved self- 


(3-TO^)A = (4>|,(r)$a:(r))-|-($]/(r)$j^(r))-|-(4't(r)4'(r)), 

(3.7) 

- ^ = (4>l(r)4-^(r)) - ($);(r)$y(r)), (3.8) 

A^ = (3.9) 

where m is magnitude of the SC condensate in SC phase, 
and we have redefined ry and g to absorb a factor of i. 
Af in Eq. (3.8) is the mean-field nematic order parameter 
defined as the anisotropy between CDWs along x and y 
directions. In other words, in our theory the nematic or¬ 
der is a vestigial order of CDW^®, breaking C 4 rotational 
symmetry but preserving lattice translational symmetry 
(see Sec. VHA for discussion of other possible origins of 
nematic order). We numerically solve the above equa¬ 
tions by computing the integrals with the given lattice 
regularization. 


B. Monte Carlo 

Classical Monte Carlo methods are capable of measur¬ 
ing standard equilibrium thermodynamic estimators of 
our model, such as energy, magnetization, CDW and SC 
structure factors, and various correlation lengths. Our 
simulations are performed on finite-size lattices and in¬ 
volve a combination of local^^“^® and non-local"*’®’^® im¬ 
portance sampling techniques, as described in detail in 
Ref. 37. Non-local sampling is especially important at 
low temperatures, where both efficiency and ergodicity 
issues can become significant. Note that our non-local 
sampling involves a modified Wolff cluster update that 
is only possible when J = K, J' = 0 and Jz — 14 in 
our model. In all of the following plots, the large-A and 
corresponding Monte Carlo results adopt the same in¬ 
put parameters and can be directly compared with each 
other, except that the large-A mean-field calculations set 
J' = 0.01 J, while the Monte Carlo calculations set J' = 0 
in order to enable cluster sampling and thus avoid non- 
ergodic behaviour. Careful studies reveal that this slight 
difference in parameters does not have a significant effect 
on the structure factors shown in our plots. 

In the presence of random-field disorder {a ^ 0), 
Monte Carlo calculations of the CDW structure factor 
5$^ (k = 0) (which will be defined in Sec. IV) require av¬ 
eraging over many independent realizations of disorder, 
{hai}. Our numerical studies reveal that, as a is in¬ 
creased, the distribution of over various Realizations 
of Disorder (ROD) becomes increasingly asymmetric due 
to the fact that 5'$^ is a complicated, non-linear function 
of the disorder fields hai- As a result, the average value, 
[(‘^$x)]rod’ distribution becomes different from its 

typical value, exp [In ( 5 '$^)]j^qj 3 , where (...) and [.. 
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FIG. 1: Zero-temperature, zero-disorder phase diagram (see 
Eq. (2.2) for definitions of g and g'). For g > 0 and g' > 0 it 
is energetically expensive to form CDW order, therefore the 
ground state is SC; for g < 0 and g' < 0 the system prefers 
unidirectional CDW. The dashed and solid lines mark first- 
and second-order transitions, respectively. A bi-critical point 
is located at g' — A, g — 0. The stripe (SC-l-stripe) phase 
becomes a nematic (SC-fnematic) phase in the presence of 
disorder. The stars mark the three sets of inpnt parameters 
used in onr calculations. 


denote thermal and disorder averages, respectively. How¬ 
ever, in order to allow comparison with large-iV results 
(for which calculations of the typical value are extremely 
difficult), all of the following Monte Carlo results corre¬ 
spond to average values of the disorder distributions. The 
qualitative behaviour of the structure factors is similar if 
one instead examines the typical values. 

Note that, in cases where no disorder is present (cr = 
0), the error bars in our Monte Carlo results correspond 
to thermal averaging. In the presence of disorder, er¬ 
ror bars instead correspond to the standard deviation of 
the mean over independent ROD. Our results average 
over between 10^ and 10^ ROD. We find that both when 
we increase a and when we study temperatures near the 
structure factor peak, more ROD are required in order 
to obtain high-quality numerical results. 

Unless otherwise stated, Monte Carlo simulations are 
performed on lattices of size 32 x 32 x 8. By examining the 
behaviour of the structure factor 5'$^(k = 0) on larger 
lattices for a selective set of input parameters, we have 
determined that this size is generally sufficient to ensure 
that the data has converged within a few percent (at 
worst) of the infinite-size limit. 


IV. CDW STRUCTURE FACTOR 

Our starting point is the T = 0 zero-disorder phase 
diagram shown in Fig. 1. Three phases emerge from a 
bi-critical point: an SC, a stripe, and a coexisting SC 
and stripe phase. For finite T, finite disorder and non¬ 
zero interlayer coupling, the bi-critical point and phase 



FIG. 2: The CDW structure factor Sc>„,(k = 0) as a function 
of T with cr = 0, increasing Jz, and the parameters given 
in Eq. (4.2) in Region 1. Solid lines without indicated data 
points are the large-N saddle-point results, while the points 
with error bars (which are smaller than the data points in this 
case) are data from Monte Carlo. The lines connecting Monte 
Carlo data points are merely guides to the eyes. Short vertical 
lines mark the locations of SC transition temperatures, Tsc- 
For Jz = 0, one can still calculate the SC transition temper¬ 
ature from the Monte Carlo data, but there is no transition 
at mean-field level. In this case, the short vertical line corre¬ 
sponding to the large-A data instead marks the temperature 
at which the SC mass term becomes exponentially small. In 
the inset, we demonstrate that Tsc approaches Tmax from be¬ 
low as Jz is increased. Monte Carlo (MC) results for Tsc are 
consistently higher than the corresponding large-A mean-field 
(MF) results, while both MC and MF give the same estimates 
for Tmax within error. 


boundaries shift; see Sec. V for a detailed discussion of 
the evolution of the phase diagram. Using both large-V 
saddle-point methods and Monte Carlo simulations, we 
compute the CDW structure factor 


S^Jk = 0) 


l($i(k)<i>,(k)) 


(4.1) 


as a function of T using three sets of input parameters, 
as indicated by the stars in Fig. 1. S'$^(k = 0) and 
(k = 0) represent X-ray scattering intensities due to 
CDW at wave vectors Qx and Qy respectively. In all of 
the following plots of (k = 0) vs. T, Monte Carlo and 
large-V saddle-point results agree qualitatively, but show 
significant quantitative differences, especially for temper¬ 
atures close to those at which the CDW peak height is 
maximal. These differences can be reduced if 1/JV cor¬ 
rections are included (see Appendix B for details). 
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FIG. 3: S<s>^ (k = 0) as a function of T in Region 1 , with 
fixed Jz ^ 0 and increasing disorder, with parameters given 
in Eq. (4.2) but with g = 0.7. We show both large-(undeco¬ 
rated lines) and Monte Carlo (points with error bars) results. 
Short vertical lines indicate SC transition temperatures. 


A. Region 1 

1. Zero disorder (cr = 0) and various interlayer 
couplings, Jj. The input parameters of Eq. (2.2) are 
taken to be 

K = J, J' = O.OIJ (J' = 0 in Monte Carlo), 

14 = J^, g = 0.35J, g' = -0.033J, A = 0.033J. (4.2) 

This reproduces (approximately) the fitting parameters 
used in Ref. 35: 

A = l, w = -0.2, 3 = 0.35, g' = 0. (4.3) 

(See Appendix A for definitions of X,w,g and g'.) 

As shown in Fig. 2, the CDW structnre factor 
grows with decreasing temperature down to a non¬ 
zero Tmax, at which point it attains maximum value. 
Below Tmax, >S'$,„(k = 0) decreases until it reaches 
zero at r = 0 . <S'$^(k = 0 ) is convex on both sides 

of Tmax- As the system becomes more 3D-like, the 
prominence of the maximum is enhanced, and the SC 
transition temperature Tgc approaches Tmax (see inset 
of Fig. 2). Large-iV calculations find that T^c exceeds 
Tlnax when Jz > 5J (not shown). Monte Carlo simula¬ 
tions for Jz = 0 , 0.01 and 0.1 are performed on lattices 
of sizes 64 x 64, 32 x 32 x 8 , and 24 x 24 x 14, respectively. 

2. Disorder is present (a > 0) with fixed interlayer 

coupling. We fix Jz = O.OIJ, increase the value of g to 
g = 0.7J in order to ensure that the ground state is SC, 
and otherwise keep the input parameters the same as in 
Eq. (4.2). As shown in Fig. 3, the CDW structure factor 
as a function of temperature still has a similar shape, but 



FIG. 4: S<s,^ (k = 0) as a function of T in Region 2, using 
large-A (undecorated lines) and Monte Carlo (points with 
error bars) methods. Here we fix Jz = O.OIJ and increase the 
disorder strength. Input parameters are given in Eq. (4.4). 
Short vertical lines correspond to SC transition temperatures. 


for cr > 0 it develops a non-zero value at T = 0, which can 
be understood as a consequence of the disorder-pinning 
effect of CDW fluctuations. As tr increases, Tgc and T^ax 
get closer to each other, but T^c remains smaller than 


B. Region 2 

We change the values of g' and A so that the input pa¬ 
rameters for Eq. (2.2) become 

K = J, J' = O.OIJ (J' = 0 in Monte Carlo), 

Vz = Jz = O.OIJ, g = 0.35J, g' = 0.533J, A = 0.167J. 

(4.4) 

Although Region 2 and Region 1 both have SC as the 
zero-disorder ground state, their CDW structure factors 
behave very differently under the effect of disorder. As 
shown in Fig. 4, the feature of maximum intensity at 
Tmax is suppressed by disorder, in contrast with Fig. 3. 
Moreover, the structure factor begins to increase again 
as T approaches zero, unlike the situation in Region 1 
where the structure factor decreases monotonically as T 
is decreased below T^ax- In both Regions 1 and 2 the 
intensity at T = 0 is enhanced by disorder, which again 
is a disorder-pinning effect. 
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FIG. 5: S<s>^ (k = 0) as a function of T in Region 3, with 
fixed Jz ~ O.OIJ, increasing disorder, and parameters given 
in Eq. (4.5). Short vertical lines and dots mark SC and ne¬ 
matic transition temperatures, respectively. Tsc remains al¬ 
most constant when disorder strength a is varied, while the 
nematic transition temperature depends strongly on a. We 
show only large-results here, since Monte Carlo results are 
difficult to obtain with these input parameters in the presence 
of such strong disorder. 


C. Region 3 

The input parameters of Eq. (2.2) are taken to be 

K = J, J' = O.OIJ (J' = 0 in Monte Carlo), 

J^ = 14 = O.OIJ, g = -2.8J, g' = 0.667J, A = 0.167J. 

(4.5) 

In the absence of quenched randomness, there is a hnite- 
temperature transition in this region to a CDW (stripe) 
phase, and SC and stripe order coexist at low T . In the 
presence of quenched randomness, no long-range CDW 
order occurs, but for weak enough randomness, there re¬ 
main finite-T transitions below which nematic order and 
SC develop sequentially. 

As shown in Fig. 5, CDW correlations are greatly en¬ 
hanced below the nematic transition in the preferred di¬ 
rection, and these correlations grow monotonically to¬ 
wards T = 0. In contrast, the SC transition, which oc¬ 
curs at a lower temperature, has very little influence on 
the behaviour of the CDW structure factor. While in 
Figs. 2, 3 and 4 we found that the SC transition had a 
dramatic effect on CDW correlations, in this regime we 
find instead that nematicity plays an overwhelming role 
in determining the T-dependence of CDW correlations. 
In other words, SC and nematic transitions tend to de¬ 
crease and increase CDW correlations respectively, and 
nematicity always wins when these two factors compete. 
We should emphasize that this is not a fine-tuning effect; 
as long as there is a nematic phase with a critical temper¬ 
ature larger than the SC transition temperature, the lack 
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cr > CTc 
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0 

Isotropic ^ 


FIG. 6: Phase diagram for the Hamiltonian in Eq. (2.2) at 
T — 0 and cr > 0 (see text for definition of a^), with J' = 
Jz = O.OIJ. Dashed (sketched by hand) and solid lines mark 
first- and second-order transitions respectively. 


of a maximum in the thermal evolution of CDW struc¬ 
ture factor is generally observed for a variety of input 
parameters. 

We did not study the region where the ground state is 
purely stripe in the absence of disorder (the green region 
in Fig. 1). Due to the lack of SC, this region is probably 
less relevant to experiments. 


V. PHASE DIAGRAMS 

In this section we discuss in detail how the phase 
diagram evolves with increasing temperature and disor¬ 
der. All phase diagrams are determined by the large-A 
saddle-point method. 


A. Zero temperature 

As mentioned in Fig. 1, for zero disorder there are 
three phases in the g' — g phase diagram and a bi-critical 
point at {g',g) = (A,0). As shown in Fig 6, for small 
but non-vanishing disorder, the stripe (SC-|-stripe) phase 
is replaced by a nematic (SC-|-nematic) phase, and the 
position of the bi-critical point shifts continuously. There 
is a critical disorder strength Uc above which there is 
no nematic phase, and a first-order transition separates 
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the SC and isotropic phases. (See Appendix Cl for a 
discussion of how the phase boundaries are determined 
in those plots.) 


B. Finite temperature 

Consider the regimes of Fig. 1 with g' < A and g' > A, 
where A controls the tendency to break C 4 symmetry. 
We further limit our discussion to g' A and g' ^ A 
to avoid the complication of a drifting bi-critical point 
in the presence of increasing disorder as indicated in 
Fig. 6. We plot in Figs. 7 and 8 the T — g phase di¬ 
agram for fixed g' and increasing disorder in these two 
regimes. Appendix C 2 gives detailed information about 
how phase boundaries and multicritical points are deter¬ 
mined in these plots. 

For g’ <C A, the first-order transition between nematic 
and SC phases persists up to ctc, at which point the bi- 
critical point and the nematic phase disappear simulta¬ 
neously. As disorder is further increased beyond the 
SC phase continues to exist, but this phase gets pushed 
steadily towards larger g. Quenched disorder tends to pin 
CDW locally, which indirectly suppresses the SC order. 
However, this effect is mitigated for larger g, at which 
point CDW order is suppressed and SC is favored. 

For g' ^ A, the phase diagram has a tetra-critical 
point. At the critical disorder strength Uc, this tetra- 
critical point and the nematic phase simultaneously van¬ 
ish, similar to the situation of g' ^ A. Note that Fig. 7 
and Fig. 8 have the same Tmuiticriticai (temperatures cor¬ 
responding to the multicritical points) and CTc (see Ap¬ 
pendix C 2 for details). 


VI. NEMATIC CORRELATION LENGTH 

As shown in Sec. IV C, when nematic order exists at 
T = 0, the CDW structure factor S'$,„(k = 0) only 
reaches a maximum value at T = 0. In other words, when 
the nematic correlation length ^nem —>■ 00 we never ob¬ 
serve Tniax > 0. Here we turn to a regime where T^ax > 0 
and no nematic phase occurs, and study the ^nem- Specif¬ 
ically, we choose input parameters corresponding to Re¬ 
gion 1 in Fig. 1. As illustrated in Fig. 9, ^nem "C ^cdw for 
weak disorder, where ^cdw is the CDW correlation length. 
As the disorder strength increases, ^nem grows and even¬ 
tually becomes comparable with ^cdw at low temperature. 
However, ^nem still remains smaller than ^cdw Implica¬ 
tions of this result are discussed in Sec. VHA 2, and 
details of the calculation can be found in Appendix D. 


T 



T 





FIG. 7: Evolution of the T — g phase diagram for fixed 
g' ^ A, as a function of a (see (2.2) for definitions of g and 
g'). Dashed (sketched by hand) and solid lines mark first- and 
second-order transitions, respectively. The bi-critical point 
moves toward large g as disorder increases, whereas the zero- 
T nematic-SC transition point behaves non-monotonically. 


VII. DISCUSSION 

A. Relation to experiments 

1. X-ray scattering. The idea of calculating CDW 
structure factors from a non-linear sigma model and 
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FIG. 8: Evolution of the T — g phase diagram for fixed 
g' ^ A, as a function of a (see (2.2) for definitions of g and 
g'). Solid lines mark second-order phase transitions. The 
location of the tetra-critical point behaves non-monotonically 
as a increases, first moving towards smaller g until a = a^, 
and then moving towards larger g. 


comparing with X-ray data was initiated in Ref. 35, 
where a model similar to Eq. (2.2) in a 2D, disorder-free 
system was shown to give good quantitative agreement 
with X-ray data. However, two discrepancies remained. 
First of all, as T —)■ 0, the structure factor calculated 
using this model vanished, unlike what is seen in X-ray 



FIG. 9: CDW and nematic correlation lengths with interlayer 
coupling Jz = 0.01 and parameters given in Eq. (4.2). These 
calculations use the large-A saddle-point method. Dotted 
lines mark the SC transition temperatures. 


scattering experiments^^“^^d 7 . is, 20 - 22 ^ Secondly, Tgc 
was found to occur relatively far below the location of 
the CDW structure factor’s maximum, whereas X-ray 
experiments find that Tgc > In our 

model, in the presence of quenched disorder (cr 7 ^ 0 ), the 
CDW structure factor sustains a finite value at T = 0 
(due to a pinning effect). In addition, both the disorder 
and interlayer coupling present in our model prove to be 
effective for bringing Tgc closer to Tmax- 

2. Nematicity. Macroscopically, long nematic cor¬ 
relations are observed experimentally in the pseudogap 
regime of several different cuprate materials^’^^’^®^^*’’^^. 
These correlations could result from a nematic phase with 
infinite correlation length, or from a C' 4 -symmetric phase 
with strong nematic fluctuations and a finite correlation 
length that is long compared to ^cdw The parameter 
regime that leads to the best agreement with X-ray data 
(Region 1 in Fig. 1), however, does not host a nematic 
phase, nor does it have a considerably longer nematic cor¬ 
relation length than ^cdwi as shown in Fig. 9. Meanwhile, 
in the regime we have studied where there is a low-T ne¬ 
matic phase (Region 3 in Fig. I), the CDW structure 
factor constantly increases as T decreases, with no sign 
of turning down where the SC transition occurs. There 
are at least three possible explanations for these discrep¬ 
ancies: 

(i) In our model we have taken Jz = 14 for simplicity, 
but in reality the SC order below Tsc is three-dimensional, 
whereas the CDW order always remains essentially two- 
dimensional (in low magnetic fields). It is possible that 
for 14 S> Jz, where the 3D coupling makes the onset of 
SC order more robust and mean-held like, that a sharp 
depression of the CDW order could occur even where 
the nematic transition temperature is greater than the 
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FIG. 10: Integrated intensity from both large-A'^ mean-field 
(undecorated lines) and Monte Carlo (points with error bars) 
calculations. Dotted lines correspond to SC transitions. 


superconducting 

(ii) Another limitation of our model comes from the con¬ 
straint imposed in Eq. (2.4). This constraint is justi¬ 
fied at temperatures well below any mean-field ordering 
temperature^®. However, at higher temperatures, the 
mean squared magnitudes of both the SC and the CDW 
order surely diminish. Specifically, these local magni¬ 
tudes refer to the k-integrated correlation functions, 

( 7 . 1 ) 

1 r d®k — I- 

(7.2) 

and Eq. (2.4) implies that Isc + Icdw,x + Icdw,y = 3, inde¬ 
pendent of T. These integrated intensities are plotted in 
Eig. 10 as functions of temperature for input parameters 
from Region 1 of Eig. 1. Here we see that, local supercon¬ 
ducting and CDW orders happily coexist at temperatures 
above their ordering temperatures. The competition be¬ 
tween the two orders occurs predominantly as long-range 
correlations arise, at which point the system is forced to 
select one form of order or the other. In this sense, at 
least in our model, the competition does not primarily 
concern the amplitudes of the orders, which one might 
want to associate intuitively with a “pairing scale” in 
the case of superconductivity, but rather^^ involves com¬ 
petition at the level of the helicity moduli, i.e. the “su¬ 


perfluid stiffness” in the case of superconductivity. 

(iii) There is an implicit assumption in this discussion 
that the nematicity detected in experiment can be at¬ 
tributed to vestigial CDW order®®. However, nematic¬ 
ity may have other origins in the pseudogap regime of 
cuprates, which is beyond the scope of our model. This 
likely applies^® to the nematicity observed at somewhat 
lower doping concentration in YBCO, i.e. for <5 < 0.09. 

B. Conclusions 

One aspect of our study with far-reaching implica¬ 
tions is the remarkable degree to which the large-A^ 
mean-field results qualitatively - and in some cases even 
semi-quantitatively - reproduce the Monte Carlo data at 
N = 2. This observation was already apparent in Ref. 35, 
but has now been extended to a wider range of circum¬ 
stances. In particular, we have shown that this approach 
applies even in the presence of quenched disorder, where 
results obtained using this self-consistent mean-field the¬ 
ory combined with the replica trick reproduce the general 
trends seen in the disorder-averaged data from Monte 
Carlo simulations. Moving forward, this enables us to 
confidently use these approximate analytic approaches in 
other contexts. For example, we can now more fully ex¬ 
plore other effective field theories that could give rise to 
experimental features of the cuprates - and other highly 
correlated materials with complex behavior - related to 
fluctuating and intertwined orders. 
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Appendix A: Parameter conversion 


We provide here the conversion rules between our model (Eq. (2.2)) and the model used in Refs. 35 and 37. The 
Hamiltonian used in these previous studies is 


^ = -E 

{ij) xy 


2 6 ■ 

Q;=3 


-E 




2 6 

Vz ^ ^ Jz ^ ^ '^icx.'^jcx. 

0.-1 q :=3 


i: E -4+11: f i: nL) + ^ i: ik + n^r ++) s i: 


i Q;=3 


i \q;=3 


a—3 


(A.l) 


where ni+in 2 represents the SC order^arameter, and 77 , 3+1714 and 775 + mg represent the two CDW order parameters. 
This model imposes the constraint X]q=i that Refs. 35 and 37 do not include the effects of interlayer 
coupling and disorder, so that 14 = = h = 0 in these references. One can write the fields 4', $ 3 , and <i>y from 

Eq. (2.2) in discretized form in terms of the components riia as 


5'i = X (77ii,77i2, ■ ■ ■,ntN)’^ , 

$a;i = VSN X (77i,Ar+i,77i^Ar+2, . . . , 77i^2Ar)^ , 
^yi = VSN X (77i,2Ar+l, »^i,2Ar+2, • ■ • , ni^sN)’^ , 


u ^ 
hi = —— X 


^ hi^ 


N+2, ■ 


T-i,3N 


(A.2) 


where N is the number of components of each order parameter; A^ = 2 in Refs. 35 and 37 as well as in the Monte 
Carlo calculations within this paper. With the transformation defined in Eq. (A.2), it is then straightforward to show 
that the conversion rules between the parameters in Eqs. (2.2) and (A.l) are 


T = T/6, X = J/K, J, = J,, V, = V,, g = g-i{J/K-l), g' = 3{g' + A), w =- 6 A, 


(A.3) 


where T and T are the temperatures corresponding to Eqs. (2.2) and (A.l), respectively. 


Appendix B: Effect of 1/A correction 

1/A corrections do not change high and low temperature behaviors of the CDW structure factor 5'$,„(k = 0), but 
these corrections do have a quantitative effect on its maximum. In Fig. 11, we compare Monte Carlo results for 
5'$^(k = 0) to large-A results with and without 1/A correction for two different sets of parameters. In the first plot, 
we find a large effect of the 1/A corrections on the height of the peak: this is the region where the inverse propagator 
of $ is the smallest, and so small “self-energy” corrections can have a large effect on the correlator. In the second 
plot, we see that the 1/A correction greatly improves the agreement with Monte Carlo. 


Appendix C: Details of Section V 
1. Zero temperature 


For T = 0 and cr = 0, only the on-site terms in Eq. (2.2) remain and we are left with a spatially homogeneous 
ground state that satisfies 







A 





(C.l) 


where V is the volume of the system, and we have set J = K so that the nearest-neighbour terms in Eq. (2.2) become 
constants. The problem amounts to searching for minima of a two-variable function 


H{x,y) = g{x + y) 



yf + 


I 

A 


i.x + yf, 


(C.2) 
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FIG. 11: The effect of 1/A^ corrections in comparion with Monte Carlo data for two different sets of parameters. Both 
parameters sets have K = J, Vz = Jz = 0, J' = 0 and a = 0. The parameters for both plots correspond to Region 1 of Fig. 1. 


with A > 0 and constraints 


2 /> 0 , X + y < 3N. 


The results are (as shown in Fig. 1): 
• When g' > A, 


^: 2/)min 


H{0,0) 

H{-gN/[2{g' - A)],0) 

H{m, 0 ) 


(SC) for 5 > 0 

(stripe+SC) for 6 (A — g') < g < 0 
(stripe) for 5 < 6 (A — 5 '). 


• When g' < A, 


^ y)min 


H(3N,0) (stripe) for g < 3(A — g') 

i/( 0 , 0 ) (SC) forg>3(A-g')- 


(C.3) 


(C.4) 


(C.5) 


For r = 0, and a ^ 0 (Fig. 6 ) we must numerically solve Eqs. (3.7), (3.8) and (3.9) for the ^'-dependence of g, 
under the following conditions (where —4J — 2 — g — Ag' + g + g is the mass term of SC order parameter): 


• Nematic to SC-|-nematic: 

-4J - 2J^ - g - Ag' + g + fj, = 0, 

• SC to SC-hnematic: 

-4 J - 2J^ - g - Ag' + g + g. = 0, 
A7 = 0. 


= 0, T = 0; 

solve (3.7), (3.8) 

II 

0 

II 

0 

solve (3.7), (3.9), 


and (3.9). 

and d/dAf over both sides of (3.8) with 


2. Finite temperature 

The second-order phase transitions in Figs. 7 and 8 are obtained by numerically solving Eqs. (3.7), (3.8) and (3.9) 
under corresponding conditions: 

• Isotropic to SC: 

—4J — 2 J 2 — g — Ag' + g + g = 0, Af = 0, fh? = 0, solve (3.7) and (3.9). 

• Isotropic to stripe (isotropic to nematic if cr 7 ^ 0): 

A/" = 0, = 0, solve (3.7), (3.9) and d/dAf over both sides of (3.8) with Af = 0. 
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• Stripe to SC+stripe (nematic to SC+nematic if cr 7 ^ 0): 

—4J — 2 J 2 — 5 — 45 '+ ry +/X = 0, = 0, solve (3.7), (3.8) and (3.9). 

• SC to SC+stripe (SC+nematic if cr 7 ^ 0): 

—4J — 2Jz — <7 — 45 ' + ry + /X = 0, Af = 0, solve (3.7), (3.9) and d/dN over both sides of (3.8) with Af = 0. 


To find the multicritical points in Figs. 7 and 8 , we must impose the conditions 

- 4J - 2J,, - 5 - 4g + 7y + yi = 0, A7 = 0, = 0, (C. 6 ) 

and solve (3.7), (3.9) and d/dN over both sides of (3.8) with Af = 0. 

The reason Tmuiticriticai is independent of g and g' is that Eq. (3.9) is decoupled from Eqs. (3.7) and (3.8) at the 
multicritical point. At cr = CTc the situation is similar: at T = 0 Eq. (3.9) becomes trivial and we only need to solve 
Eqs. (3.7) and (3.8). 


Appendix D: Nematic correlation function 


We present some technical aspects of computing nematic correlation functions using the replica trick. Neglecting 
the terms related to superconductivity (yy and 'k), we can rewrite the replicated Hubbard-Stratonovich Hamiltonian 
of Eq. (3.1) in fc-space as 


i?replica= f ^ | A(k)4>t , (k) , (k) + H (k) (k) (k) 

Jd-k „ 


+ A7a(k) / $l^„(k + p)$a,x(p) - 4’I,j,(k + p)$a.y(p) 

J dp ^ 

/ Ef-^)‘J’a(k)$.(k), 

a.b V J / 


A7a(k)A7a(-k)iV 

4A 


(D.l) 


wliere fak = 10^, and 

A,B(k) = —2J{coskx + cosky) + 2J'{coskx — cosky) — 214 cosfc^ + yt. (D.2) 

Our goal is to compute (A/’o(k)A/’a(—k)) in the limit M ^ 0, where M is the total number of replicas. First, we 
integrate out 4) and obtain an effective Hamiltonian for N, which satisfies 




= / 


= / 


axp /?F7replica ^ 

exp I - ■ • •. , ^M,x)^ + ( 4 >k,..., $ky) Ba/k) ($l.y, ..., 4>M.y)^ I 


xexpi-/ 3E / ■^a(k) 4>l_„(k + p)$a.:.(p) - 4>),_j^(k + p)$a,y(p) 

I a •'d.kdp 


Na(k)Na{-k)N 

4A 


(D.3) 


where 


' T 


1 ... 1 

1 ... 1 


Act, Bcr(k) = A, B(k) ■ ImxM 


(D.4) 
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WVW' + 


k+p 


FIG. 12: The lowest order diagrams corresponding to Eq. (D.3) (both are order 1/N). The wavy lines correspond to the bare 
propagator of Af, and the solid lines correspond to the bare propagator of • Each vertex is order 1 and the loop is of order 
N (coming from the sum of all the "F’s). 


We now expand the last exponential term in Eq. (D.3) to quadratic order in JV (which is equivalent to keeping the 
diagrams up to order 1/A^, as shown in Fig. 12), and arrive at 

exp I - ^I.x(k + P)$a.^(p) - $I,y(k + p)$a,!/(p) 

= f ^I.:r:(k + p)$a.x(p) - $I,;;(k + p)$a.y(p) 

J d]<.dp ^ 




A/;(k)A4(k') 


^ J dk.dpd\^'dp' 

^I.:r(k + P)$a.x(p) - «'I,j,(k + p)$a,y(p) ’ „ (k' + p')$M (p') “ (k' + p')®6,!/(p') (D-5) 

The term linear in /3 will vanish due to the cancellation between $ 3 , and $j,, and we will drop the constant 1. The 
remaining /3^ term gives 


j3H,sW]=Yl AA,(k)U(k)A4(-k), 

a,b 


(D. 6 ) 


(D.7) 


u(k) = (^U^M(k + p))($^$a,x(p)) + («>l,j^«>b,y(k + p))($];_y$a.y(p)), 

where we have included from Eq. (D.3), and dropped the factor of TV. The nematic correlation function 

is 

(AA,(k)AA,(-k)) = ^ ^(A4(k)A4(-k)) = ^(sr' + ... + sl^) (D.8) 


where {sa} are eigenvalues of 2tab (the factor 2 is due to the fact that the $'s are complex). Our strategy is to first 
compute ($t 2 .<i)(,_a;(k)) in Eq. (D.7) by diagonalizing ^^(k), which gives 




T 


M’A{\l) M A{k)-2a^M/T' 
We then diagonalize tab and obtain (A/’a(k)A/’Q(—k)) according to Eq. (D. 8 ), which yields 
(AA„(k)AA,(-k))(M^0) 


(D.9) 


2AT Jip *1^ 


1 I 2 (T^ 

T(k+p)A(p) + T 


T2(k+p)T(p) ^ A(k+p)A2(p) 


_ loT r 
T2 Jc 


dp A2(k+p)A2(p) 


+ A^B 


^ 2AT Sdp *1^ A{k+p)A(p) 


20-2 

T 


A2(k+p)A(p) ^ A(k+p)A2(p) 


+ A-H-B 


I) 


At k = 0, Eq. (D.IO) reduces to 


(AA„(0)A/;(0))(M^0) = 


2AT Jdp ^ 






A2(p) w B^{p) ^ TA3(p) ^ TB3(p) 


) 


( 4A 


dp \ T2A-‘(p) + T2B4 (p 


4cr'‘ 


1 

2AT 


fdp ( A2( 


A2(p) + B2(p) TA3(p) TB^(p 


4(7^ _|_ 4(7^ 


(D.IO) 


(D.ll) 
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where A and B are defined in Eq. (D.2). A nematic transition occurs when the denominator vanishes such that 

1 1 4 ct 2 40-2 


1 


2AT 


V^"(P) ' S2(p) ' TA3(p) ' Ti?3(p)) 


(D.12) 


which is consistent with mean-field saddle-point equation (3.8) after taking a derivative with respect to M and setting 
A/" = 0 in (3.8). 
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